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Abstract 

The  purpose  of  this  paper  is  to  describe  a  method  for  simulation  of  recently  introduced 
fluid  stochastic  Petri  nets.  Since  such  nets  result  in  rather  complex  set  of  partial  differential 
equations,  numerical  solution  becomes  a  formidable  task.  Because  of  a  mixed,  discrete  and 
continuous  state  space,  simulative  solution  also  poses  some  interesting  challenges,  which  are 
addressed  in  the  paper. 
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1  Introduction 


Stochastic  Petri  nets  provide  a  convenient  and  concise  method  of  describing  discrete-event 
dynamic  systems  [1,  4,  6,  12,  15].  One  of  the  difficulties  encountered  while  using  stochastic 
Petri  nets  is  that  the  underlying  reachability  graph  tends  to  be  very  large  in  practical  prob¬ 
lems.  Drawing  a  parallel  with  fluid  flow  approximations  in  performance  analysis  of  queueing 
systems  [3,  7,  11],  SPNs  have  been  extended  to  include  fluid  (or  continuous)  places  where 
fluid  can  be  used  to  approximate  a  large  number  of  discrete  tokens.  Armed  with  such  Fluid 
Stochastic  Petri  nets  (FSPNs),  we  can  also  model  physical  systems  that  contain  continuous 
fluid-like  quantities  which  are  controlled  by  discrete  logic. 

Fluid  stochastic  Petri  nets  were  introduced  by  Trivedi  and  Kulkarni  in  [14].  The  original 
model  was  considerably  enhanced  in  [9].  The  purpose  of  this  paper  is  twofold:  first  we 
further  extend  the  formalism  in  [9]  to  make  it  more  useful  and,  second,  we  explore  the  use 
of  simulation  as  a  solution  method  for  FSPNs.  The  extensions  to  the  FSPN  formalism  we 
propose  include: 

•  Fluid  impulses  associated  with  both  immediate  and  timed  transition  firings. 

•  Guards  on  immediate  transitions,  dependent  on  fluid  levels,  not  just  on  the  discrete 
marking. 

In  the  process  of  using  simulation  as  a  solution  method,  several  innovations  are  needed. 
In  this  direction,  the  contributions  of  this  paper  include: 

•  Generation  of  random  deviates  based  on  a  non-homogeneous  Poisson  process. 

•  Interleaving  of  ODE  solution  for  fluid  places  with  simulation  of  discrete  events  in  the 
FSPN. 

•  Definition  of  restrictions  under  which  one  can  integrate  the  change  of  fluid  levels  using 
built-in  closed-from  results,  such  as  decoupled  behavior  and  special  classes  of  functions 
for  the  fluid  rates. 

After  introducing  the  FSPN  model  in  the  next  section,  we  describe  the  method  of  sim¬ 
ulation  for  the  most  general  case  in  Section  3.  The  simplest  case  of  uncoupled  behavior  of 
different  fluid  places  is  taken  up  in  Section  4  while  two  other  cases  are  discussed  in  Sections 
5  and  6,  respectively.  Examples  are  provided  in  Section  7. 
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2  Fluid  stochastic  Petri  nets 


Sets  are  denoted  by  upper  case  calligraphic  letters.  In  particular,  M ,  IZ,  and  TZo  indicate 
the  natural,  real,  and  nonnegative  numbers,  respectively. 

For  simplicity,  we  only  address  exponentially  distributed  firing  times.  Generally  dis¬ 
tributed  firing  times  are  clearly  useful,  and,  in  connection  with  discrete-event  simulation, 
might  not  add  much  complexity  to  the  solution.  However,  the  interruption  policies  (what 
happens  to  the  remaining  firing  time  of  transitions  when  one  of  them  fires)  can  require  very 
complex  definitions  in  full  generality  [5,  13].  This  is  not  the  case  with  the  exponential  distri¬ 
bution,  due  to  its  memoryless  property.  If  we  assume  that  the  firing  rate  of  each  transition, 
if  marking-dependent,  is  reevaluated  in  each  marking  where  it  is  enabled,  the  time  elapsed 
since  the  transition  first  became  enabled  does  not  affect  the  future  evolution  of  the  FSPN. 

A  fluid  stochastic  Petri  net  (FSPN)  is  a  tuple  (V^ ,T^ ,'T^ ,«-■,  f  i Qi  )? 

where: 

•  =  {pi, . . .  ,P|^D|}  and  =  {qi, . . . ,  q\vc\}  are  two  disjoint  and  finite  sets  of  places. 

Let  V  =  U  V^.  A  (discrete)  place  p  e  h  drawn  with  a  single  circle  and  can 
contain  a  number  of  tokens  nip  G  Af.  A  (continuous)  place  q  G  is  drawn  with 
two  concentric  circles  and  can  contain  a  level  of  fluid  G  TZo-  The  marking,  or 
state,  of  the  FSPN  is  given  by  a  pair  of  vectors  describing  the  contents  of  each  place, 
(m,x)  G  <5  =  X  We  call  S  the  “potential  state  space”,  as  opposed  to 

the  “actual  state  space”  S  C  S,  the  set  of  markings  actually  reachable  during  the 
evolution  of  the  FSPN.  The  marking  (m,x)  evolves  in  time,  which  we  indicate  by  r, 
so,  formally,  we  can  think  of  it  as  a  stochastic  process  {(m(r),x(r)),r  >  0}. 

•  =  {ti, . . .  ,f|7-T|}  and  =  {ui, . . . ,  W|r^|}  are  two  disjoint  and  finite  sets  of  transi¬ 
tions.  Let  T  =  7^  U  T^.  A  (timed)  transition  t  G  7^  is  drawn  as  a  rectangle  and  has 
an  exponentially  distributed  firing  time.  An  (immediate)  transition  f  G  is  drawn 
as  a  thin  bar  and  has  a  constant  zero  firing  time. 

•  a  :  X  T)  U  (T  X  V^))  X  5  ->  and  a  :  x  T)  U  (T  x  V^))  x  <S  7^o 

describe  the  marking-dependent  cardinality  (for  discrete  places)  or  the  fluid  impulse 
(for  continuous  places)  of  the  input  and  output  arcs  connecting  transitions  and  places. 
We  use  the  same  symbol  for  both,  and  we  draw  them  an  thin  arcs  with  an  arrowhead 
on  their  destination,  since  the  type  of  place  eliminates  any  possibility  of  confusion. 
The  function  describing  a  is  written  on  the  arc,  the  default  is  the  constant  one. 
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•  /  :  {{V^  X  T)  U  (T  X  V^))  X  5  -¥  'IZo  describes  the  marking-dependent  fluid  rate  of 
the  input  and  output  arcs  connecting  transitions  and  continuous  places.  These  fluid 
arcs  are  drawn  with  a  thick  line,  and  an  arrowhead  on  their  destination.  Also  in  this 
case  the  function  is  written  on  the  arc  and  the  default  is  the  constant  one. 

•  g  :  T  x<S— )-{0, 1}  describes  the  marking-dependent  guard  of  each  transition. 

•  A  :  X  S  TZo  and  w  :  X  >  TZo  describe  the  marking-dependent  firing  rates 

and  weights  of  each  transition. 

•  b:V^xAf\^^\  TZqU  {00}  describe  the  fluid  bounds  on  each  continuous  place.  This 

bound  has  no  effect  when  it  is  set  to  infinity.  Note  that  b  depends  only  on  the  discrete 
part  of  the  state  space,  not  on  <5,  to  avoid  the  possibility  of  circular  definitions. 

•  (m°,  x“)  G  is  the  initial  marking.  Graphically,  it  is  represented  by  writing  the  value 
of  irip,  or  x°,  inside  the  corresponding  place.  A  missing  value  indicates  zero.  For 
discrete  places,  it  is  also  common  to  draw  nip  tokens  inside  the  place,  if  this  number 
is  small. 

The  meaning  of  these  quantities  is  given  by  the  enabling  and  firing  rules.  We  say  that  a 
transition  f  G  T  has  concession  in  marking  (m,x)  iff 

Vp  G  ap,t(m,x)  <  nip  and  pt(ni,x)=l. 

If  any  immediate  transition  has  concession  in  (m,  x),  it  is  said  to  be  enabled  and  the  marking 
is  said  to  be  vanishing.  Otherwise,  the  marking  is  said  to  be  tangible  and  any  timed  transition 
with  concession  is  enabled  in  it.  In  other  words,  a  timed  transition  is  not  enabled  in  a 
vanishing  marking  even  if  it  has  concession. 

Some  definitions  of  SPNs  allow  one  to  disable  a  transition  t  with  concession  in  a  marking 
by  specifying  a  zero  rate  or  weight  for  it,  or  by  introducing  inhibitor  arcs,  drawn  with  a  circle 
instead  of  an  arrowhead.  Since  we  can  represent  these  behaviors  by  an  appropriate  definition 
of  input  arc  cardinalities  or  of  guards,  we  assume,  without  loss  of  generality,  that  rates  and 
weights  are  positive  for  an  enabled  transition.  Inhibitor  arcs  can  then  be  considered  merely 
as  a  shorthand^. 

§If,  in  (m,x),  an  inhibitor  arc  from  p  G  {q  G  V^)  to  t  E  T  has  cardinality  c  E  N  [c  E  t  is 
disabled  if  c  >  (c  >  Xg).  The  same  behavior  can  be  modeled  by  ensuring  that  the  guard  gt  evaluates  to 
0  in  (m,  x). 
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Let  £^(ni,x)  denote  the  set  of  enabled  transitions  in  marking  (m,x).  Enabled  transitions 
change  the  marking  in  two  ways.  First,  a  transition  t  ^'T  enabled  in  marking  (ni,x)  can  fire 
after  a  random  amount  of  time  having  distribution  ~  Expo(A((m,x)),  and  yield  a  (possibly) 
new  marking  (m',  x').  We  then  write  (m,x)-4(m',  x'),  where 

\/peV^,  =  mp  + at,p(m,x)  -  ap,t(m,x) 

Vq  e  x'  =  min{6,(m'),max{0,X5  +  ai,,(m,x)  -  a5,t(m,x)}}. 

Second,  fluid  flows  continuously  through  the  arcs  /  of  enabled  transitions  connected  to 
continuous  places.  The  potential  rate  of  change  of  fluid  level  for  the  continuous  place  q  G 
in  marking  (m,  x)  is 


(yf*(m,x)=  Yl  /M(m,x)-/5,t(m,x). 

t£S{xn,x) 

However,  the  fluid  level  can  never  become  negative  or  exceed  the  bound  65,  (m),  so  the  (actual) 
rate  of  change  over  time,  r,  while  in  marking  (m,x),  is 


5g{m,x)  =  ^ 


0  if  bq{m)  =  0  or 

if  x,  =  bq{m)  and  (5^®‘(m,x)  >  0  or 
if  Xg  =  0  and  (^P°*(m,x)  <  0 
x)  otherwise 


(1) 


The  stochastic  evolution  of  the  ESPN  in  a  tangible  marking  is  governed  by  a  race  [2]: 
the  timed  transition  t  with  the  shortest  firing  time  is  the  one  chosen  to  Are  next,  unless 
it  becomes  disabled  by  some  fluid  levels  reaching  particular  values  that  cause  t  to  become 
disabled  prior  to  its  firing.  In  a  vanishing  marking,  instead,  the  weights  are  used  to  decide 
which  transition  should  Are:  an  immediate  transition  u  enabled  in  marking  (m,  x)  fires  with 
probability 

w„(m,x)  ^2) 

Y  t«„'(m,x) 


3  General  case 

The  ESPN  definition  we  just  gave  is  very  powerful,  but  it  allows  one  to  describe  models 
whose  solution  can  be  quite  difficult,  even  with  discrete-event  simulation.  Indeed,  it  can  be 
used  to  define  FSPNs  whose  behavior  is  “unstable,”  as  in  the  FSPNs  of  Fig.  1.  In  the  model 
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Figure  1:  FSPNs  exhibiting  unstable  behaviors. 

on  the  left,  immediate  transitions  Ui  and  U2  alternatively  put  and  remove  a  unit  impulse 
instantaneously.  With  few  exceptions  [8],  such  a  behavior  has  been  considered  a  modeling 
error  in  the  literature  on  discrete-state  models.  The  instability  of  the  model  in  the  middle 
is  instead  exclusive  to  models  with  a  states  having  a  continuous  component,  such  as  our 
FSPNs.  When  x®  =  0,  timed  transition  ti  is  enabled  and  timed  transition  ^2  is  disabled. 
However,  as  soon  as  the  fluid  arc  starts  adding  fluid  to  q,  the  situation  is  reversed,  ti  becomes 
disabled,  while  ^2  becomes  enabled  and  starts  emptying  q.  It  could  be  argued  that,  in  such  a 
situation,  q  will  always  be  empty,  but  any  model  where  an  infinite  number  of  events  occurs  in 
a  finite  time  (e.g.,  transitions  ti  and  ^2  become  enabled  an  infinite  number  of  times)  cannot 
be  managed  by  conventional  discrete-event  simulation  techniques.  Hence,  we  will  consider 
such  behaviors  illegal. 

The  model  on  the  right  could  be  also  considered  unstable  if  F2  >  Fi.  Both  ti  and  ^2  are 
always  enabled,  hence  there  is  a  continuous  flow  into  q  at  rate  Fi  due  to  ti.  However,  the 
outgoing  flow  due  to  ^2  cannot  be  F2.  Our  definition  simply  states  that  5q  is  0  in  this  case, 
implying  that  the  outgoing  flow  is  effectively  reduced  to  be  Fi,  instead  of  Fj,  or,  in  other 
words,  the  arc  from  q  to  ^2  can  be  thought  to  have  effect  only  a  fraction  F1/F2  of  the  time. 
This  type  of  behavior,  however,  can  be  easily  managed  by  considering  all  the  flows  incident 
to  a  continuous  place,  so  we  do  not  regard  it  as  a  true  instability. 

We  now  describe  how  to  study  a  model  with  no  unstable  behaviors.  Assume  that  we 
have  just  entered  tangible  marking  (m,x).  If  there  is  any  enabled  transition,  each  continues 
component  x,  might  vary  in  a  very  general  way  over  time.  Applying  Eq.  1  to  each  q  €  we 
obtain  a  system  of  ordinary  differential  equations  subject  to  the  initial  condition  x(0)  =  x. 
We  can  then  consider  two  cases: 

•  In  the  simpler  case,  the  cardinality  of  the  arcs  connected  to  discrete  places  and  the 
guards  do  not  depend  on  x.  Even  so,  the  firing  times  behave  as  a  nonhomogeneous 
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Poisson  process  (NHPP)  whose  rate  depends  on  the  continuous  marking,  and  so  some 
care  is  required  in  sampling  the  firing  instants.  We  assume  that  the  firing  rate  of 
each  transition  t  can  be  bounded  from  above  by  A*(m),  given  our  knowledge  of  its 
dependence  on  the  fluid  marking.  That  is,  when  the  discrete  marking  is  m,  the  rate  of 
t  satisfies  At(m,x)  <  A^  (m),  for  any  value  of  x  that  might  be  reached  in  conjunction 
with  m.  We  can  therefore  sample  from  the  NHPP  using  the  technique  of  “thinning” 
[10],  where  we  sample  “potential  firing  instants”  in  accordance  with  a  homogeneous 
Poisson  arrival  process  with  rate 

A*(m)=  E 

From  this  process,  we  can  define  a  sequence  of  increasing  time  instants 
Starting  from  i  —  we  accept’’  that  is,  we  declare  that  a  firing  occurred  at  time 
Ti,  with  probability  A(m,  x(ri))/A*(m),  where 

A(m,x(ri))=  Y.  A<(m,x(Ti)). 

In  other  words,  we  use  the  actual  firing  rates  at  time  as  a  weight,  to  determine 
whether  the  event  corresponds  to  a  true  firing  or  not.  This  requires  us  to  solve  for 
the  value  of  x(ri),  by  integrating  the  system  of  ordinary  differential  equations.  If  n  is 
accepted,  we  stop.  Otherwise,  we  integrate  until  r2,  compute  x(r2),  and  decide  whether 
to  accept  r2  or  not,  and  so  on.  Eventually,  this  process  stops  at  some  step  i,  giving  us 
the  actual  firing  time  =  Ti, 

For  example.  Fig.  2  illustrates  the  case  where  four  transitions  are  enbaled  in  (ni,x),  ti, 
^2,  is,  and  The  sequence  of  numbered  arrows  shows  the  random  deviates  generated, 
in  order.  First,  we  generate  Ti  according  to  the  distribution  Expo(A*(m)).  Then  we 
generate  a  random  deviate  ~  Unif(0,  A*(m)).  In  the  figure,  this  happens  to  fall  in 
the  interval  corresponding  to  the  “do  not  accept”  case.  Thus,  we  need  to  generate 
another  potential  firing  time  by  sampling  the  distribution  Expo(A  (ni))  again  and 
summing  the  sampled  value  to  rj,  obtaining  T2-  We  also  need  another  random  deviate 
~  Unif(0,  A*(m)),  which  also,  in  the  figure,  happens  to  cause  a  rejection.  Finally, 
we  generate  a  third  potential  firing  time,  When  we  sample  ~  Unif(0,  A*(m)),  we 
obtain  a  value  falling  in  the  interval  corresponding  to  ^2,  hence  we  schedule  the  firing 
of  t2  at  time  It  is  then  apparent  that  the  expected  number  or  random  deviates  that 
need  to  be  generated  is  larger  when  the  bounds  At(m,x)  for  the  enabled  transitions 
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Figure  2:  Sampling  the  NHPP  process  underlying  a  FSPN. 

are  less  tight.  On  the  other  hand,  if  the  rates  of  the  enabled  transitions  are  a  function 
of  X,  but  Z]<€£(in,x) is  a  known  constant  independent  of  x,  only  two  deviates 
are  needed:  the  first  one  to  decide  Tj  and  the  second  one  to  decide  which  transition  to 
fire. 

•  If,  instead,  the  set  of  enabled  transitions  can  change  as  x  evolves,  we  also  need  to 
consider  an  “enabling  event”  at  the  time  r®  when  the  first  change  in  £’(m,  x)  occurs. 
The  method  to  compute  r®  depends  on  the  nature  of  the  dependencies.  In  principle, 
we  should  know  the  value  of  x(t)  over  the  entire  horizon  r  G  [0,r-^].  This  can  still 
be  accomplished  through  integration.  After  (during)  integration  we  need  to  find  that 
value  of  r®  that  first  satisfies  the  given  condition  on  the  fluid  levels.  If  there  is  no 
minimum  value  r®  €  [0,  for  which  the  set  of  enabled  transitions  changes,  the  next 
event  to  schedule  is  the  firing  at  time  Otherwise,  we  must  schedule  an  “enabling 
event”  at  time  r®. 

In  either  case,  if  the  firing  rates  of  timed  transitions  are  not  dependent  on  fluid  levels,  the 
generation  of  next  firing  times  is  considerably  simplified  because  the  machinery  of  NHPP- 
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based  generation  of  random  deviates  is  avoided. 

The  processing  of  the  scheduled  event  causes  a  change  of  marking,  from  (m,  x)  to  (m  ,  x  ), 
where  =  m  if  the  event  was  of  the  enabling  type.  Then,  in  marking  (m^,  x  ),  a  finite 
sequence  of  immediate  firings  might  take  place,  just  as  in  ordinary,  non-fluid,  SPNs,  until 
the  next  tangible  marking  (m",x")  is  reached.  Thanks  to  the  memoryless  property  of  the 
exponential  distribution,  the  evolution  of  the  FSPN  from  this  point  on  is  analogous  to 
the  evolution  from  the  initial  marking,  that  is,  we  do  not  need  to  be  concerned  about 
the  “remaining  firing  times”  of  transitions  that  were  already  enabled  prior  to  reaching  this 
marking. 


4  Uncoupled  behavior 

The  general  behavior  just  described  requires  us  to  solve  a  system  of  ordinary  differential 
equations  at  each  step  of  the  simulation.  This  computation  can  be  quite  costly.  A  restriction 
on  the  type  of  dependency  allows  us  to  uncouple  the  system,  resulting  in  a  set  of  ordinary 
differential  equations  which  can  be  solved  independently.  This  requires  that  the  fluid  rates 
incident  on  hence  (^^(m,x),  depend  only  on  (m,Xg),  not  on  the  fluid  levels  in  the  other 
continuous  places: 

V(m,x),(m,x')  eS,  =  x'  ^g(m,x)  =  ^^(m,x'). 

As  in  the  general  case,  we  can  still  distinguish  whether  the  set  of  enabled  transitions  can 
be  affected  by  x  or  not,  and  the  NHPP  random  variate  generation  must  be  used  only  if  their 
firing  rates  depend  on  x. 


5  Predefined  classes  of  behaviors 


For  particular  cases  of  uncoupled  dependencies,  we  can  even  have  a  built-in  closed  form 
solution,  which  will  avoid  the  need  for  numerical  integration  altogether.  One  such  case  is 
when,  in  a  given  naarking  (m,x), 

=  A(m)x,(T)  +  B(m),  A(m)  /  0 


whose  solution  is 


x,(t) 


B(m) 

>l(m) 


+  x,(0)  + 


B{m)\ 

A{m)j 
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assuming  that  x,  remains  between  0  and  ^>g(ni)  during  [0,r].  This  answers  the  question  of 
how  much  the  fluid  level  in  a  place  will  change  during  the  firing  time  r  of  a  timed  transition. 
Inversely,  the  time  r,  when  place  q  reaches  a  certain  fluid  level  threshold  Lg  is  given  by 


In 


if  this  quantity  is  positive  (if  it  is  negative,  we  can  simply  define  r,  =  oo,  that  is,  the 
threshold  Lg  cannot  be  reached  in  this  marking). 

If  the  set  of  enabled  transitions  can  only  change  when  some  place  q  reaches  a  threshold 
level  Lg,  then  we  can  simply  define  the  time  r®  of  the  next  enabling  event  as 

T®  =  min{rg}. 

When  A(m)  =  0,  the  solution  is  much  simpler, 

— =  B(m)  ^  Xg(r)  =  Xg(0)  +  .B(m)r, 

dr 

again  assuming  that  Xg  remains  between  0  and  feg(m)  during  [0,r].  The  time  Tg  when  place 
q  reaches  the  threshold  Lg  is  then 

_  Lg  -  Xg(0) 

^  “  B{m)  ’ 

if  this  quantity  is  positive,  infinity  otherwise. 

6  Piecewise  constant  behavior 

Complete  dependency  on  the  marking  (m,  x)  is  desirable  in  principle,  but  the  complication 
it  entails  is  often  excessive  and  its  full  power  unneeded.  A  simpler  type  of  dependency  is 
obtained  by  enforcing  a  discretization  on  the  behavior  of  the  FSPN  with  respect  to  the 
continuous  component  x.  This  can  be  accomplished  by  defining  a  set  of  boolean  threshold- 
type  conditions  C  =  {(ri  ©i  L), . . . ,  (r|£|  ©|£|  /|£|)},  where  is  a  continuous  place, 

Qk  G  {<5  ^5  =1  ^5  >»  is  a  comparison  operator,  and  Ik  :  — >  P  U  {oo}  is  a  threshold 

dependent  on  the  discrete  marking  only.  Hence,  given  a  marking  (m,x),  we  can  define 


/!(m) 


Xg(0)+ 


H(m) 


A(m) 
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the  “discretized”  marking  (m,  i)  obtained  through  £,  where  i  G  {0, and  ife  —  1  iff 
Xr^  Qk  4(ni). 

If  we  force  a  (for  discrete  places  only),  /,  g,  and  A  to  be  defined  on  the  discretized 
marking  (m,i),  rather  than  on  (m,x),  then  the  behavior  of  the  FSPN  is  constant  until  the 
first  threshold  is  encountered,  or  until  a  firing  occurs. 

Hence,  we  can  carry  on  a  traditional  discrete-event  simulation,  where  the  types  of  events 
that  need  to  be  scheduled  in  the  event  queue  are  either  transition  firings  or  the  hitting  of  a 
threshold. 

Fortunately,  there  is  no  need  to  place  the  same  restriction  on  the  fluid  impulses  (a  for 
continuous  places)  or  w,  since  the  impulses  and  the  weights  are  always  evaluated  only  at  a 
specific  instant  in  time.  Applying  the  restriction  to  these  quantities  as  well  would  prevent  us 
from  modeling  useful  behaviors,  such  as  emptying  a  continuous  place,  or  choosing  between 
two  immediate  transitions  with  probability  proportional  to  the  level  in  two  continuous  places. 


7  Examples 

We  illustrate  the  power  of  the  formalism  with  a  few  examples. 

7.1  A  queue  with  impatient  customers  and  breakdowns 

Consider  a  queue  with  a  server  subject  to  breakdowns  and  repairs.  The  customers  arrive 
with  a  constant  rate,  and  queue  in  an  unbounded  waiting  room.  They  are  served  in  first- 
come-first-serve  order,  but,  once  their  service  starts,  they  can  become  impatient  and  leave 
before  completion,  see  Fig.  3.  Unlike  other  system  with  impatient  customers,  the  amount  of 
time  a  customer  has  been  in  the  queue  before  his  service  begins  does  not  affect  his  decision 
to  leave.  The  arcs  from  Serving  to  Busy  and  from  Waiting  to  Idle  are  used  to  count  time 
into  the  two  places,  hence  they  have  fluid  rate  one.  The  arcs  from  Busy  and  Idle  to  Serving 
have  impulse  xsusy  and  Xidie  defined  on  them,  respectively.  Hence,  they  are  “flushing”  arcs, 
they  have  the  effect  of  emptying  the  two  places  immediately  after  the  firing  of  Serving. 

The  guard  of  immediate  transition  Leave  specifies  when  the  customer  at  the  head  of  the 
queue  decides  to  leave.  Various  policies  can  be  easily  modeled: 

•  The  total  amount  of  time  from  the  moment  service  began  exceeds  a  certain  threshold 
MAX.  Then,  we  could  define  the  guard  gLeave  to  be  the  boolean  expression  {xBu$y  + 
y^Idle  =  MAX). 
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Figure  3;  The  FSPN  of  a  queue  with  impatient  customers  and  breakdowns. 

This  policy  is  representative  of  situations  where,  once  the  server  begins  operating  on 
a  customer,  the  operation  must  complete  within  a  certain  time,  to  avoid  spoilage,  etc. 

•  The  total  amount  of  time  a  customer  has  not  received  any  service  from  the  moment 
service  began  exceeds  a  certain  threshold  MAX.  Then,  gLeave  =  {^idie  =  MAX). 

This  could  represent  a  similar  situation,  where,  spoilage  occurs  only  when  the  customer 
is  not  being  served. 

•  A  customer  has  waited  for  an  uninterrupted  period  of  time  MAX  without  receiv¬ 
ing  any  service.  Then,  gieave  =  {'X-uie  —  MAX),  after  adding  an  an  impulse  arc 
asusj/,Repair(ni,x)  =  Xsusy,  SO  that  place  Busy  becomes  empty  after  each  repair. 

This  could  represent  a  situation,  where,  in  addition  to  occurring  only  when  the  cus¬ 
tomer  is  not  being  served,  any  spoilage  immediately  disappears  as  soon  as  service 
resumes. 

•  A  customer  has  spent  more  time  waiting  for  the  server  to  be  operational  than  receiving 
service,  from  the  moment  service  began.  Then,  gieave  =  {^idie  >  '^-Busy)- 
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A  measure  of  interest  is  the  fraction  of  jobs  that  decide  to  leave: 

number  of  firings  of  Leave  up  to  time  r 
number  of  firings  of  Arrive  up  to  time  r 

computed  over  a  finite  horizon,  or  in  the  limit  for  r  oo. 

7.2  A  dual-tank  processing  facility 

Consider  a  processing  plant  where,  during  normal  operation,  a  liquid  enters  a  main  tank. 
One,  from  an  external  source  with  rate  7i„,  and  is  used  by  a  processing  station,  with  a 
(potential)  rate  'font  >  7m ,  (see  Fig.  4). 

However,  the  processing  station  is  subject  to  breakdowns  during  which  it  cannot  process 
the  liquid.  Interrupting  the  flow  from  the  external  source  of  liquid  into  the  main  tank  is  an 
expensive  operation,  hence,  a  second  additional  tank.  Two,  is  present.  When  the  processing 
station  is  down,  the  liquid  is  automatically  routed  to  tank  Two,  which  has  a  maximum 
capacity  bxwo-  Only  when  the  second  tank  is  full,  the  flow  from  the  external  source  is  shut 
down.  After  a  repair,  the  processing  can  resume  and  the  liquid  is  routed  again  from  the 
external  source,  which  is  restarted  if  it  had  been  shut  down,  into  tank  One.  In  addition,  any 
liquid  in  tank  Two  is  pumped  into  tank  One,  with  rate  712.  If  7m  +  712  >  7oiit?  th®  level 
in  tank  One  will  increase  while  the  processing  station  is  working  to  catch  up  after  a  repair. 
Since  tank  One  has  a  maximum  capacity,  bone-,  the  flow  from  tank  Two  to  tank  One,  rather 
than  the  flow  from  the  external  source,  is  slowed  down  when  this  limit  is  reached.  The  guard 
(in  the  FSPN  of  Fig.  5)  gxjer  =  (xone  <  bone)  accomplishes  this. 

The  main  reason  for  having  two  tanks,  instead  of  a  single  large  one,  is  efficiency.  As  the 
liquid  needs  to  be  maintained  at  a  given  temperature,  tank  One  is  constantly  heated,  while 
tank  Two  is  heated  only  when  it  contains  liquid,  i.e.,  during  a  breakdown.  Indeed,  the  two 
measures  we  are  interested  in  computing  are: 

number  of  firing  of  Start  up  to  time  r 

■  ^ 

T 

the  frequency  at  which  the  external  source  needs  to  go  through  a  start-stop  cycle,  and 
probability  that  tank  Two  is  not  empty  at  time  r, 
again,  either  for  a  finite  r  or  in  the  limit  for  r  — >■  00. 
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Figure  4:  A  dual-tank  processing  facility. 


Figure  5:  The  FSPN  of  the  dual-tank  processing  facility. 

8  Conclusion 

In  this  paper  we  have  extended  the  power  of  recently  introduced  fluid  stochastic  Petri  nets. 
Since  equations  characterizing  the  evolution  of  such  FSPNs  are  a  coupled  system  of  partial 
differential  equations,  the  generation  and  solution  of  these  equations  can  become  intractable 
except  for  small  or  very  well  structured  FSPNs.  Hence,  discrete-event  simulation  becomes 
an  important  avenue  for  the  solution  of  FSPNs.  Because  of  mixed,  discrete  and  continuous 
state  space  with  heavy  interactions  between  them,  simulation  also  poses  some  challenges. 
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Some  of  the  challenges  are  addressed  in  the  paper.  Actual  implementation  (currently  in 
progress)  of  an  FSPN  simulator  will  undoubtedly  reveal  other  problems. 


References 

[1]  M.  Ajmone  Marsan,  G.  Balbo,  and  G.  Conte,  Performance  Models  of  Multiprocessor 
Systems,  MIT  Press,  Cambridge,  MA,  1986. 

[2]  M.  Ajmone  Marsan,  G.  Balbo,  A.  Bobbio,  G.  Chiola,  G.  Conte,  and  A.  Cumani,  “The 
effect  of  execution  policies  on  the  semantics  and  analyis  of  Stochastic  Petri  Nets,”  IEEE 
Trans.  Softw.  Eng.,  vol.  15,  pp.  832-846,  July  1989. 

[3]  D.  Anick,  D.  Mitra  and  M.  Sondhi,  “Stochastic  Theory  of  Data-Handling  Systems,”  The 
Bell  System  Technical  Journal,  Vol.  61,  No.  8,  pp.  1871-1894,  Oct.  1982. 

[4]  C.  G.  Cassandras,  Discrete  Event  Systems:  Modeling  and  Performance  Analysis,  Aksen 
Associates,  Holmwood,  IL,  1993. 

[5]  G.  Ciardo,  “Discrete-time  Markovian  stochastic  Petri  nets,”  in  Numerical  Solution  of 
Markov  Chains  ’95  (W.  J.  Stewart,  ed.),  (Raleigh,  NC),  pp.  339-358,  Jan.  1995. 

[6]  G.  Ciardo,  J.  K.  Muppala,  and  K.  S.  Trivedi,  “Analyzing  concurrent  and  fault-tolerant 
software  using  stochastic  Petri  nets,”  J.  Par.  and  Distr.  Comp.,  15(3):255-269,  July 
1992. 

[7]  A.  I.  Elwalid  and  D.  Mitra,  “Statistical  Multiplexing  with  Loss  Priorities  in  Rate-Based 
Congestion  Control  of  High-Speed  Networks,”  IEEE  Transaction  on  Communications, 
Vol.  42,  No.  11,  pp.  2989-3002,  November  1994. 

[8]  W.  K.  Grassmann  and  Y.  Wang,  “Immediate  events  in  Markov  chains,”  in  Numerical 
Solution  of  Markov  Chains  ’95  (W.  J.  Stewart,  ed.),  (Raleigh,  NC),  pp.  163-176,  Jan. 
1995. 

[9]  G.  Horton,  V.  Kulkarni,  D.  Nicol,  and  K.  Trivedi,  “Fluid  stochastic  Petri  nets:  Theory, 
application,  and  solution,”  ICASE  Report  96-5,  Institute  for  Computer  Applications  in 
Science  and  Engineering,  Hampton,  VA,  1996. 


14 


[10]  P.A.W.  Lewis  and  G.S.  Shedler,  “Simulation  of  Nonhomogeneous  Poisson  Processes  by 
Thinning,”  Naval  Research  Logistics  Quarterly,  Vol.  26,  pp.  403-414,  1979. 

[11]  D.  Mitra,  “Stochastic  Theory  of  Fluid  Models  of  Multiple  Failure-Susceptible  Producers 
and  Consumers  Coupled  by  a  Buffer,”  Advances  in  Applied  Probability,  Vol.  20,  pp.  646- 
676,  1988. 

[12]  T.  Robertazzi,  Computer  Networks  and  Systems:  Queueing  Theory  and  Performance 
Evaluation,  Springer- Verlag,  1990. 

[13]  M.  Telek,  A.  Bobbio,  and  A.  Puliafito,  “Steady  state  solution  of  MRSPN  with  mixed 
preemption  policies,”  in  Proc.  IEEE  International  Computer  Performance  and  Depend¬ 
ability  Symposium  (IPDS’96),  (Urbana-Cliampaign,  IL,  USA),  pp.  106-115,  IEEE  Comp. 
Soc.  Press,  Sept.  1996. 

[14]  K.  S.  Trivedi  and  V.  C.  Kulkarni,  “FSPNs:  fluid  stochastic  Petri  nets,”  in  Proc.  l^th 
Int.  Conf.  on  Applications  and  Theory  of  Petri  Nets,  (Chicago,  IL),  pp.  24-31,  June 
1993. 

[15]  N.  Viswanadham  and  Y.  Narahari,  Performance  Modeling  of  Automated  Manufacturing 
Systems,  Prentice-Hall,  Englewood  Cliffs,  NJ,  1992. 


15 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporUng  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  inforniation,  including  suggestions  for  reducingthis  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


1.  AGENCY  USE  ONLYf/.eai/e  blank)  2.  REPORT  DATE 

_ _  May  1997 


4.  TITLE  AND  SUBTITLE 

DISCRETE-EVENT  SIMULATION  OF  FLUID 
STOCHASTIC  PETRI  NETS 


6.  AUTHOR(S) 

Gianfranco  Ciardo 
David  Nicol 
Kishor  S.  Trivedi 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science  and  Engineering 
Mail  Stop  403,  NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


3.  REPORT  TYPE  AND  DATES  COVERED 

Contractor  Report 


5.  FUNDING  NUMBERS 

C  NASl-19480 
WU  505-90-52-01 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  97-24 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23681-0001 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

NASA  CR- 201688 
ICASE  Report  No,  97-24 


11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 
Final  Report 

(Accepted)  PNPM  ^97,  St.  Malo,  France,  June  1997 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Unclassified-Unlimited 


12b.  DISTRIBUTION  CODE 


Subject  Category  60,61 


13.  ABSTRACT  (Maximum  200  words) 

The  purpose  of  this  paper  is  to  describe  a  method  for  simulation  of  recently  introduced  fluid  stochastic  Petri 
nets.  Since  such  nets  result  in  rather  complex  set  of  partial  differential  equations,  numerical  solution  becomes 
a  formidable  task.  Because  of  a  mixed,  discrete  and  continuous  state  space,  simulative  solution  also  poses  some 
interesting  challenges,  which  are  addressed  in  the  paper. 


14.  SUBJECT  TERMS  15.  NUMBER  OF  PAGES 

discrete-event  simulation,  continuous  system  simulation,  Petri  Nets,  Stochastic  Petri  17 

N^ets 

16.  PRICE  CODE 

_ _ _ _ _  A03 


17,  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION 

OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT  OF  ABSTRACT 

U  nclassified  Unclassified 


NSN  7540-01-280-5500 


Standard  Form  298(Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


